# This is an example how the trackanalysis is performed for one experiment

# For each experiment, three types of input files are necessary

# A) Information about the environment / about the pores of the chip (same for all examples)
# CFW-TL-RL-zstack-4x4-merged-with-FITC-2_Mai-modified_Pores.xls

# B) Imaris segmentation
# CFW-TL-RL-zstack-4x4-merged-with-FITC-2_Mai-modified_Hyphae_semi-manual_3Filaments.wrl

# C) Information about filaments
# CFW-TL-RL-zstack-4x4-merged-with-FITC-2_Mai-modified_Hyphae_semi-manual_3Filaments.xls

# NOTE 1: In my setting, the raw data are one folder below the current working directory ("../")
# NOTE 2: In the filenames of the raw data, spaces (" ") were replaced by a single dash ("-") 

# 1. This data must be preprocessed now:

# Using the cpp code ReadWRL (which must be compiled beforehand), we transform the .wrl file into a .csv file:
# On command line:
# $ ReadWRL/cmake-build-debug/readWRL ../CFW-TL-RL-zstack-4x4-merged-with-FITC-2_Mai-modified_Hyphae_semi-manual_3Filaments.wrl processed_files/segmentation_3_hyphae.csv

# Using in2csv, we transform the .xls files into csv files extracting the necessary sheets:
# On command line:
# Firstly, the pores
# $ in2csv ../CFW-TL-RL-zstack-4x4-merged-with-FITC-2_Mai-modified_Pores.xls --sheet 'Position' | tail -n +2 >> processed_files/pores.csv
# Secondly, the information about filaments  
# $ in2csv ../CFW-TL-RL-zstack-4x4-merged-with-FITC-2_Mai-modified_Hyphae_semi-manual_3Filaments.xls --sheet 'Dendrite Branch Level' | tail -n +2 >> processed_files/filament_information_3.csv

# Everything saved to the folder 'processed_files' - this is already done! You can start from here!

# Start analysis
# Set the working directory to the folder where this script is located (individual on each computer)
#setwd("~/Code_Analysis_Confinement_Ratio")
# Set folder of processed files
folder = "processed_files/"

# Load data
set = read.csv(paste(folder, "segmentation_3_hyphae.csv", sep=""),sep = ' ',header = FALSE)
fila = read.csv(paste(folder, "filament_information_3.csv", sep=""),header = TRUE)

### 2. Approximation for the membrane 

# An algorithm is applied which creates a grid approximating the membrane by the position of the (2 closest) pores
gridsize = 5 # other values possible - the smaller the better the accuracy, but higher computation costs

# load pores
pores = read.csv(paste(folder, "pores.csv", sep=""))

# assisting function
roundTo = function(x,y){ ##round to multiple of given y
  if((y - x %% y) <= x %% y) { x + (y - x %% y)}
  else { x - (x %% y)}
}

#get boundaries of the grid
minx = roundTo(min(as.numeric(pores$Position.X),set$V3)-gridsize,gridsize)
maxx = roundTo(max(as.numeric(pores$Position.X),set$V3)+gridsize,gridsize)
miny = roundTo(min(as.numeric(pores$Position.Y),set$V4)-gridsize,gridsize)
maxy = roundTo(max(as.numeric(pores$Position.Y),set$V4)+gridsize,gridsize)

## The grid is the same for all three examples and can be directly loaded below!
regGrid = data.frame(x=pores$Position.X,y=pores$Position.Y,z=pores$Position.Z)  ## the grid always contains also the pores.
#compute the z level for each grid point
for(xL in seq(minx-gridsize,maxx+gridsize,gridsize)){
  print(paste0(xL))
  for(yL in seq(miny-gridsize,maxy+gridsize,gridsize)){
    ##compute xy-distance of xL,yL to all the pores
    dist = sqrt((pores$Position.X-xL)**2+(pores$Position.Y-yL)**2)
    p2 = pores
    p2$dist = dist
    p2$weights = max(dist)-dist
    p2 = p2[with(p2, order(-weights)), ]
    p2 = p2[1:2,] #take the closest 2 points
    wm = weighted.mean(x = p2$Position.Z,w = p2$weights)# make a weighted mean of their z level as appoximation
    #        print(wm)
    regGrid = rbind(regGrid,c(xL,yL,wm))# collect output
  }
}
# save the calculated grid for the membrane
#write.csv(x = regGrid,file = paste(folder, "grid.csv", sep=""))

## Load grid
regGrid = read.csv(paste(folder, "grid.csv", sep=""),header = TRUE)

### 3. Calculate new z-level (based on the membrane)
colnames(set) = c("FilamentID","ignore","x","y","z","ignore2")
set$xR = roundTo(set$x,gridsize) # the warnings can be ignored
set$yR = roundTo(set$y,gridsize) # the warnings can be ignored
set$pore_level_aprox = Inf

set$xR[set$xR < 0.1 & set$xR > -0.1] <- 0
set$yR[set$yR < 0.1 & set$yR > -0.1] <- 0

#subset(regGrid, regGrid$x==50 & regGrid$y==50)
for(xL in seq(minx-gridsize,maxx+gridsize,gridsize)){
  #print(paste0(xL))
  for(yL in seq(miny-gridsize,maxy+gridsize,gridsize)){
    #get the corrected membene level
    z_approx = subset(regGrid,regGrid$x==xL & regGrid$y==yL)$z
    set$pore_level_aprox[set$xR==xL & set$yR==yL] = z_approx
    #print(paste0(z_approx))
    if(z_approx > 1000) {
      print(paste0(xL, yL))
    }
  }
}
set$dist2membraneAprox = set$z-set$pore_level_aprox # distance for all hyphae to approximated membrane
#now, set contains the new 'dist2membraneAprox'. this is basically the new z coordinate we use now!

### 4. Add the filament information to the data set

set$HyphaeID = -1
set$depth = -1
set$level = -1
# Add the information in the filament file to the data set 
for(i in unique(set$FilamentID)){ 
  if(dim(subset(fila,fila$ID==i))[1]>0){
    set$HyphaeID[set$FilamentID==i] = subset(fila,fila$ID==i)$FilamentID
    set$depth[set$FilamentID==i] = subset(fila,fila$ID==i)$Depth
    set$level[set$FilamentID==i] = subset(fila,fila$ID==i)$Level 
  }
}

# Print some information about the data set
print(paste0("The data set contains ", length(unique(set$HyphaeID))," Hyphae and ", length(unique(set$FilamentID)), " filaments"))
# 3 hyphae
# 351 filaments

# We can have a look at the result, different hyphae are colored
library("plot3Drgl", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.4")
plot3D::scatter3D(x = set$x,y=set$y,z = set$dist2membraneAprox,colvar = set$HyphaeID)
plotrgl()

### 5. Find piercing branches
# Now, piercing filaments and their belonging hyphae are identified
# Then, all side branches which do not constitute to the connection between conidium and the membrane piercing part are removed to apply a skeletonization later

filaIDs = vector()
hyphaeIDs = vector()
for(i in unique(set$FilamentID)){
  a = subset(set,set$FilamentID==i)
  if( max(a$dist2membraneAprox) > 0 & min(a$dist2membraneAprox) < 0){
    #print(paste0("Filament ",i," is going through the membrane"))
    filaIDs = c(filaIDs, i)
    hyphaeIDs = c(hyphaeIDs, unique(a$HyphaeID))
  }
}
print(paste0("There are ", length(unique(filaIDs)), " piercing filaments of ", length(unique(hyphaeIDs)), " hypha"))
# 4 piercing filament of 1 hypha

# Since we have more than one piercing filament, we have to loop over all
resA = data.frame(x=0,y=0,z=0,N=0,ID=0)
for (fId in 1:length(filaIDs)) {
  print(paste0("Process branch ", fId))
  # Start analysis of a single filament and find its branch recursively to the conidium
  
  chosen_branch = filaIDs[fId]
  set2 = subset(set,set$FilamentID==chosen_branch)
  
  ### 6. Begin skeletonization
  a = set2
  resBranch = a
  # recursively collect all segments from here back to the conidium.
  nowDepth = unique(a$depth)
  print(paste0("Start: nowDepth:",nowDepth))
  # We only need 60 µm of the branch to reconstruct for the analysis, therefore a higher depth level than 0 is sufficient
  depth_necessary = max(0, nowDepth - 4) # this value was chosen for each data set individually
  while(nowDepth != depth_necessary){
    # find all candidates
    #print(paste0("nowDepth:", nowDepth))
    b = subset(set,set$depth==(nowDepth-1))  # look in the full data for all filaments with a - 1 depth
    # find the one which is closer
    min = 100000;
    minID = NA
    minDepth = NA
    # now we have to find the filament which is closest to our current one
    for(j in 1:nrow(a)){
      b$dist = sqrt((a[1,]$x-b$x)**2+(a[1,]$y-b$y)**2+(a[1,]$z-b$z)**2)
      # find min in b$dist and get the corresponding Sub_ID and save it.
      smallest = b[b$dist==min(b$dist),]
      if(smallest$dist < min){
        min = smallest$dist
        minID = smallest$FilamentID # save the ID of the closest filament
        minDepth = smallest$depth
      }
    }
    
    a = subset(set,set$FilamentID==minID) # now this is the current one
    resBranch = rbind(resBranch,a)
    nowDepth = unique(a$depth) # set new before re entering the loop
  }
  
  # Now, the branch is reconstructed
  # plot3D::scatter3D(x = resBranch$x,y=resBranch$y,z = resBranch$dist2membraneAprox,colvar = resBranch$FilamentID)
  # plotrgl()
  
  # We can apply the sliding window function to get a skeletonization
  with = 50 # these parameters were found sufficient
  step = 50
  resBranch = resBranch[with(resBranch, order(depth)), ]
  library("zoo", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.4") 
  xM = rollapply(resBranch$x, width = with, by = step, FUN = mean, na.rm = TRUE, align = "left")
  yM = rollapply(resBranch$y, width = with, by = step, FUN = mean, na.rm = TRUE, align = "left")
  zM = rollapply(resBranch$dist2membraneAprox, width = with, by = step, FUN = mean, na.rm = TRUE, align = "left")
  out = data.frame(x=xM,y=yM,z=zM)
  
  # plot3D::scatter3D(x = out$x,y=out$y,z = out$z,colvar = out$z)
  # plotrgl()
  
  # In the next steps, the quality of the skeletonization will be increased
  # In the end, the idea is to have a point each 1 µm
  res = data.frame(x=0,y=0,z=0,N=0,ID=0)
  dat = out
  lp = dat[1,]
  length = 0;
  dist = 1 
  counter = 1
  for(i in 2:length(dat$x)){
    #get the new distance the hyphae has grown
    preal = dat[i,]
    addDist = sqrt( (preal$x-lp$x)**2+(preal$y-lp$y)**2+(preal$z-lp$z)**2 )
    length = length+addDist
    #print(paste0("last point is :",round(lp$x,2),"|",round(lp$y,2),"|",round(lp$z,2),")"))
    #print(paste0("At point ",i," the track has length",round(length,3), "(i:",round(preal$x,2),"|",round(preal$y,2),"|",round(preal$z,2),")"))
    if(length > dist){
      #print(paste0("inserting new points!"))
      # need to add a new (or mulitple ones) to the new track
      # how many point to i need to add?
      nr = floor(length/dist)
      for(j in 1:nr){
        # and now add the point
        sx = (preal$x-lp$x)/length # shift along x
        sy = (preal$y-lp$y)/length
        sz = (preal$z-lp$z)/length
        # and now add the new shifted point
        newPoint_x =lp$x + j*sx 
        newPoint_y =lp$y + j*sy 
        newPoint_z =lp$z + j*sz
        counter = counter + 1
        res = rbind(res,c(newPoint_x,newPoint_y,newPoint_z,counter))
        # now the new point was set,  need to update the rest
        if(j==nr){
          lp = data.frame(x=newPoint_x,y=newPoint_y,z=newPoint_z)
          length = length%%dist
        }
      }
    }
  }
  res = res[-1,]
  # Now we only take the 60 µm
  res = subset(res, res$z <= 0) # start at the membrane
  #print(paste0("Length of branch starting from membrane: ", length(res[,1]))) # Check the size of the branch so far (if smaller than 60, the depth level above has to be changed)
  length_branch = 60
  res = subset(res, res$N > max(res$N)-length_branch) # remove everything larger than 60 µm
  res$ID = fId
  
  #plot3D::scatter3D(x = res$x,y=res$y,z = res$z,colvar = res$z)
  #plotrgl()
  
  resA = rbind(resA, res)
  
}

resA = resA[-1,]

### 7. Track Analysis is applied to skeletonized hypha
# This package had been implement Marco Blickensdorf and Sören Totzauer in 2015 as 'projekt modul'
# The package is an implementation of the cell track analysis from zeinab (please read https://journals.plos.org/plosone/article/comments?id=10.1371/journal.pone.0080808)
# It takes tracks (time ordered points in 2d or 3d) and computes some measures for it: asphericity, confinement ratio, displacement ratio and outreach ratio (abbreviated CAOD)
install.packages("TrackAnalysis/CelltrackStats2.tar.gz", repos = NULL, type = "source")
library("CelltrackStats2", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.4")
# Also necessary
library("ggplot2", lib.loc="/usr/lib/R/site-library")
resA$t = resA$N
# This init function needs information to process the data
# Type "3d" in the first question and can just press return for remaining 
init()
# The package only operates on a data frame calles "track_data"
track_data = resA
#track_data$ID = 1 # the ID paramater must be set. If you want to run multiple cell tracks at once, put them all into track_data but give them different IDs
build_track_data()
# An output folder is generated, where the analysis files are located (onle confinement ratio is needed)

### 8. Look at the final plot for the confinement ratio 

# Insert output folder here
output_folder = "output_3_hyphae"

res = data.frame(Var1=0,Var2=0,value=0,ID=0,m=0)
for(i in list.files(paste0(output_folder),full.names = FALSE)){
  if(grepl("confinement", i, fixed=TRUE)){
    dat = read.csv(paste0(output_folder, "/",i))
    dat = dat[,-1]
    dat = melt(as.matrix(dat))
    length =max(dat$Var1)
    dat$Var2 = rep(1:length, each = length, times = 1)
    id = as.numeric(strsplit(i, "_")[[1]][3])
    dat$ID = id
    dat$m = "C"
    res = rbind(res,dat)
  }
}
res = res[-1,]

dat = res
# package: dplyr
plot = dat %>% group_by(Var1,Var2,m) %>% summarise(avg = mean(value)) %>% data.frame()
minval = min(plot$avg)
ggplot(subset(plot,plot$m=="C"),aes(Var1,Var2))+geom_raster(aes(fill = avg))+
  scale_fill_gradientn(colours = rainbow(200)[130:1],limits=c(min(0.7,minval),1))+
  coord_fixed()+labs(x="hyphal length [µm]",y="hyphal length [µm]",fill="Confinement Ratio")